switch experiment

    case 'sizedep'

        p.taus        = Mu - 1;           % size-dependent subsidy
    
        Xnew          = fsolve('findequilibrium', [D; Y*1.05; N*1.1], options, w, z, p, 'planner', 'new');

    case 'efficient'
    
        p.taus        = 0; 
        Xnew          = fsolve('findequilibrium', [D*1.2; Y*1.6; N*1.2], options, w, z, p, 'planner', 'new');

end

% New Steady State Allocations

Dnew          = Xnew(1); 
Ynew          = Xnew(2); 
Nnew          = Xnew(3); 

[~, qnew, ~, Onew, Znew, Wnew, Lpnew, Knew, Bnew, Unew, Uqnew, Cnew, Munew] = findequilibrium(Xnew, w, z, p, 'planner', 'new');

Gnew          = (Znew/Dnew/(1 + p.taus) - 1)/Munew;
Qnew          = 1/(1 - p.beta*(1 - p.varphi))*Ynew/Nnew*Gnew; 
Lnew          = Lpnew + p.F*p.varphi*Nnew; 
Rnew          = 1/p.beta - 1 + p.delta; 

ynew          = qnew*Ynew; 

% Now start doing transition dynamics

global cfun fspace; 

npts      = 25; 

fspace    = fundefn('cheb', npts, 0.9, 1.5); 

Ngrid     = funnode(fspace); 

Zgrid     = zeros(npts, 1); 
Ogrid     = zeros(npts, 1); 
Mugrid    = zeros(npts, 1); 
Dgrid     = zeros(npts, 1);
Pgrid     = zeros(npts, 1);

options.Display = 'none';

for i = 1 : npts

    Dp      = fsolve('findequilibrium', Dp, options, w, z, p, 'planner', 'old', Ngrid(i));

    [~, ~, ~, Ogrid(i), Zgrid(i), ~, ~, ~, ~, ~, ~, ~, Mugrid(i), Pgrid(i)]  = findequilibrium(Dp, w, z, p, 'planner', 'old', Ngrid(i));
    
    Dgrid(i) = Dp; 
    
end
    
Ggrid  = (Zgrid./Dgrid/(1 + p.taus) - 1)./Mugrid;               % flow gains from varieties

cfun   = funfitxy(fspace, Ngrid, [Zgrid, Ggrid, Mugrid]);

xpar   = [p.beta; p.varphi; p.delta; p.psi; p.F; p.nu; p.alpha; p.theta; p.phi];                       save xpar xpar; 

x1_ss  = [K;    N];                                                                                    save x1_ss x1_ss; 
x2_ss  = [Ynew; Cnew; Bnew; Lpnew; Knew; Nnew; Wnew; Rnew; Qnew; Znew; Gnew; Onew; Munew];             save x2_ss x2_ss; 

dynare transition noclearall;

% Cleanup

filename = 'transition';

eval(['delete ' filename '*.m']);
eval(['delete ' filename '*.mat']);
eval(['delete ' filename '*.swp']);
eval(['delete ' filename '*.log']);
eval(['delete ' filename '*.eps']);
eval(['delete ' filename '*.pdf']);
eval(['delete ' filename '*.fig']);
eval(['delete ' filename '*.jnl']);

% Populate first period values from original steady state

Yt(1)  = Y; 
Ct(1)  = C; 
Bt(1)  = B; 
Lpt(1) = Lp; 
Kt(1)  = K; 
Nt(1)  = N; 
Wt(1)  = W; 
Rt(1)  = R; 
Zt(1)  = Z; 
Mut(1) = Mu; 

Lt     = [L; Lpt(2:end) + p.F*(Nt(2:end) - (1 - p.varphi)*Nt(1 : end-1))]; 
Xt     = [p.delta*K; Kt(2:end) - (1 - p.delta)*Kt(1 : end - 1)]; 

T    = length(Ct) - 1; 

Wold = 1/(1 - beta)*(log(C)  - psi/(1 + nu)*L^(1 + nu)); 

Wnew = (beta.^(0 : 1 : T - 1))*(log(Ct(2 : end)) - psi/(1 + nu)*Lt(2:end).^(1 + nu)) + beta^T/(1 - beta)*(log(Ct(end)') - psi/(1 + nu)*Lt(end)'.^(1 + nu));

dW   = (exp((1 - beta)*(Wnew - Wold)) - 1)*100;

clc

switch experiment

    case 'sizedep'

        Experiment = 'Size-dependent Subsidy';

    case 'efficient'

        Experiment = 'Efficient Allocation';

end

fprintf('<strong> Table 4: Implications of Alternative Policies, Benchmark Model </strong> \n'); 

fprintf('\n');

fprintf(' Steady state comparisons in percent - aggregate markup at %1.2f \n', Mu);  

fprintf('\n');

fprintf([' Experiment --> ' Experiment, ' \n']);  

fprintf('\n');

disp(table(round((Ynew/Y   - 1)*100, 1), round((Cnew/C   - 1)*100, 1), round((Lnew/L   - 1)*100, 1), round((Nnew/N   - 1)*100, 1), round((Knew/K   - 1)*100, 1), round((Znew/Z   - 1)*100, 1), round(dW, 2), ...
    'VariableNames', {'Y','C','L','N','K','Z','Welfare,%'}))
